package cn.fxc;

import org.apache.commons.math3.distribution.FDistribution;
import org.apache.commons.math3.distribution.TDistribution;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;

import java.util.*;

public class MultipleLinearRegressionComputer {
    //因变量数组
    private double[] y;
    //自变量矩阵
    private double[][] x;
    private String[] xName;

    //系数
    private double[] beta;
    //系数标准误差
    private double[] parametersStandardErrors;
    //(回归)自由度
    private int degreesOfFreedom;
    //T值
    private double[] tValue;
    //P值
    private double[] pValue;
    //方差膨胀因子
    private double[]VIF;

    //西格玛，S
    private double sigma;
    //R平方
    private double rSquared;
    //R平方调整
    private double adjustedRSquared;
    //R平方预测
    private double predRSquared;
    //方差分析表
    private List<Map<String,Object>>ANOVAMapList;
    //显著性参考值(画图用)
    private double criticalValue;
    //自变量按t值的帕累托排序(画图用)
    private List<Parameter> paretoParameterList;

    public boolean doMultipleLinearRegression(){
        try {
            if (y==null)
                throw new RuntimeException("y is null");
            if (x==null)
                throw new RuntimeException("x is null");
            for (double[] xi:x){
                if (y.length!=xi.length){
                    throw new RuntimeException("the length of data is not same");
                }
            }

            OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
            double[][] xt = MatrixUtils.createRealMatrix(x).transpose().getData();
            regression.newSampleData(y, xt);

            beta = regression.estimateRegressionParameters();
//        for (double d:beta)
//            System.out.println("beta: "+d);

            parametersStandardErrors = regression.estimateRegressionParametersStandardErrors();
//        for (double d:estimateRegressionParametersStandardErrors)
//            System.out.println("estimateRegressionParametersStandardErrors: "+d);
            degreesOfFreedom = y.length-x.length-1;
            TDistribution tDistribution = new TDistribution(degreesOfFreedom);
            tValue = new double[beta.length];
            pValue = new double[beta.length];
            for (int i=0;i<beta.length;i++){
                tValue[i] = beta[i]/parametersStandardErrors[i];
    //            System.out.println("T Value["+i+"]: "+tValue[i]);
                pValue[i] = 2* (1- tDistribution.cumulativeProbability(Math.abs(tValue[i])));
    //            System.out.println("P Value["+i+"]: "+pValue[i]);
            }

            double SST = regression.calculateTotalSumOfSquares();
            double SSE = regression.calculateResidualSumOfSquares();

            double SSR = SST-SSE;
            double adjMS = SSR/x.length;
//        double RSq = (Math.abs(SST) < 1e-10)?(1-Double.MAX_VALUE):(1.0- (SSE/SST));
//        System.out.println("RSq: "+RSq);
            double dfResidual = y.length - x.length -1;
            rSquared = regression.calculateRSquared();
//        System.out.println("rSquared: "+rSquared);
            double[]rSquaredi = new double[x.length];
            VIF = new double[x.length];
            double[]SSRxi= new double[x.length];
            double[]adjSSRxi = new double[x.length];
            double[]fxi = new double[x.length];
            double[]pxi = new double[x.length];
            double adjMSe = SSE/dfResidual;
            double dfR = x.length;
            double fValueR = adjMS/adjMSe;
            FDistribution fDistribution = new FDistribution(x.length,y.length-x.length-1);
            double pValueR = 1.0 - fDistribution.cumulativeProbability(fValueR);
            for (int i=0;i<x.length;i++){
                double[] newY = x[i];
                double[][] newX = new double[x.length-1][];
                for (int t=0,s=0;s<x.length;s++){
                    if (s!=i){
                        newX[t]= x[s];
                        t++;
                    }
                }
                double[][] newXt = MatrixUtils.createRealMatrix(newX).transpose().getData();
                OLSMultipleLinearRegression regressionXi = new OLSMultipleLinearRegression();
                regressionXi.newSampleData(newY, newXt);
                rSquaredi[i]= regressionXi.calculateRSquared();
    //            System.out.println("rSquaredi["+i+"]:"+ rSquaredi[i]);
                VIF[i] = 1/(1-rSquaredi[i]);
    //            System.out.println("VIF["+(i+1)+"]:"+ VIF[i]);

                OLSMultipleLinearRegression regressionYi = new OLSMultipleLinearRegression();
                regressionYi.newSampleData(y,newXt);


                double SSTxi = regressionYi.calculateTotalSumOfSquares();
                double SSExi = regressionYi.calculateResidualSumOfSquares();
                SSRxi[i] = SSTxi-SSExi;
                adjSSRxi[i] = SSR-SSRxi[i];
    //            System.out.println("adjSSRxi["+i+"]"+adjSSRxi[i]);
                fxi[i]=adjSSRxi[i]/1.0/adjMSe;
                pxi[i]=1.0-fDistribution.cumulativeProbability(fxi[i]);
            }
//        System.out.println();
//
//
//        for (double d:residuals)
//            System.out.println("residuals: "+d);
//        double[][] parametersVariance = regression.estimateRegressionParametersVariance();
//        for (double[] pArray:parametersVariance) {
//            for (double d : pArray) {
//                System.out.println("parametersVariance: " + d);
//            }
//            System.out.println();
//        }
//
//        double regressandVariance = regression.estimateRegressandVariance();
//        System.out.println("regressandVariance: "+regressandVariance);

            sigma = regression.estimateRegressionStandardError();
//        System.out.println("S: "+sigma);
            rSquared = regression.calculateRSquared();
//        System.out.println("Rsq: "+rSquared);
            adjustedRSquared = regression.calculateAdjustedRSquared();
//        System.out.println("Rsq-调整: "+adjustedRSquared);
            double[][] hat = regression.calculateHat().getData();
            double press = 0;
            for(int i=0;i<y.length;i++){
                press +=Math.pow(regression.estimateResiduals()[i]/(1-hat[i][i]),2);
            }
            predRSquared = 1.0 - press/SST;
//        System.out.println("Rsq-预测: "+predRSquared);

            ANOVAMapList = new ArrayList<>();
            Map<String, Object> regressonMap = new HashMap<>();
            regressonMap.put("name", "回归");
            regressonMap.put("DOF", dfR);
//        System.out.println("回归自由度(dfM): "+dfR);
            regressonMap.put("adjustedSS", SSR);
//        System.out.println("SSR(回归adjSS): "+SSR);
            regressonMap.put("adjustedMS", adjMS);
//        System.out.println("回归adjMS: "+adjMS);
            regressonMap.put("FValue", fValueR);
//        System.out.println("回归FValue: "+fValueR);
            regressonMap.put("PValue", pValueR);
//        System.out.println("回归PValue: "+pValueR);
//        System.out.println();
            ANOVAMapList.add(regressonMap);
            for (int i=0;i<x.length;i++){
    //            System.out.println("第"+(i+1)+"个变量的自由度: 1.0");
    //            System.out.println("第"+(i+1)+"个变量的adjSS: "+adjSSRxi[i]);
    //            System.out.println("第"+(i+1)+"个变量的adjMS: "+adjSSRxi[i]/1.0);
    //            System.out.println("第"+(i+1)+"个变量的FValue: "+fxi[i]);
    //            System.out.println("第"+(i+1)+"个变量的PValue: "+pxi[i]);
                Map<String, Object> xMap = new HashMap<>();
                if (xName!=null)
                    xMap.put("name", xName[i]);
                else {
                    xMap.put("name", String.valueOf(i));
                }
                xMap.put("DOF", 1.0);
                xMap.put("adjustedSS", adjSSRxi[i]);
                xMap.put("adjustedMS", adjSSRxi[i]/1.0);
                xMap.put("FValue", fxi[i]);
                xMap.put("PValue", pxi[i]);
                ANOVAMapList.add(xMap);
            }
            Map<String, Object> residualMap = new HashMap<>();
            residualMap.put("name", "残差");
            residualMap.put("DOF", dfResidual);
            residualMap.put("adjustedSS", SSE);
            residualMap.put("adjustedMS", adjMSe);
            ANOVAMapList.add(residualMap);
//        System.out.println("残差自由度(dfResidual): "+dfResidual);
//        System.out.println("SSE(残差adjSS): "+SSE);
//        System.out.println("残差adjMS: "+adjMSe);

            TDistribution tDistributionStand = new TDistribution(y.length-x.length-1);
            criticalValue  = tDistributionStand.inverseCumulativeProbability(1 - 0.05 / 2);
//        System.out.println("criticalValue : "+criticalValue );
            paretoParameterList = new ArrayList<>();
            for (int i=1;i<=x.length;i++){
                paretoParameterList.add(new Parameter(String.valueOf(i),tValue[i]));
            }
            paretoParameterList.sort(Comparator.reverseOrder());
        } catch (RuntimeException e) {
             e.printStackTrace();
             return false;
        }

        return true;
    }

    public double[] getBeta() {
        return beta;
    }

    public double[] getParametersStandardErrors() {
        return parametersStandardErrors;
    }

    public int getDegreesOfFreedom() {
        return degreesOfFreedom;
    }

    public double[] gettValue() {
        return tValue;
    }

    public double[] getpValue() {
        return pValue;
    }

    public double[] getVIF() {
        return VIF;
    }

    public double getSigma() {
        return sigma;
    }

    public double getrSquared() {
        return rSquared;
    }

    public double getAdjustedRSquared() {
        return adjustedRSquared;
    }

    public double getPredRSquared() {
        return predRSquared;
    }

    public List<Map<String, Object>> getANOVAMapList() {
        return ANOVAMapList;
    }

    public double getCriticalValue() {
        return criticalValue;
    }

    public List<Parameter> getParetoParameterList() {
        return paretoParameterList;
    }

    public double[] getY() {
        return y;
    }

    public void setY(double[] y) {
        this.y = y;
    }

    public double[][] getX() {
        return x;
    }

    public void setX(double[][] x) {
        this.x = x;
    }

    public String[] getxName() {
        return xName;
    }

    public void setxName(String[] xName) {
        this.xName = xName;
    }


    public static void main(String[] args) {

        double[] y = new double[]{53,58,52,40,60,58,77,64,47,59,70,78,59,53,48,57,46,52,47,55,66,45,53,78,41,56,75,45,46,63,45,76,61,79,55,64,54,46,72};
        double[][] x = new double[2][];
        x[0] = new double[]{1,1,0,0,0,1,1,1,0,1,1,1,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,1,1,0,0,1};
        x[1] = new double[]{179,167,171,164,169,169,175,158,162,167,163,161,170,174,173,168,163,156,173,159,177,178,170,155,161,165,165,158,157,170,155,168,171,169,174,158,176,179,176};
        String[] name = new String[]{"性别","身高"};
        MultipleLinearRegressionComputer multipleLinearRegressionComputer =new MultipleLinearRegressionComputer();
        multipleLinearRegressionComputer.setY(y);
        multipleLinearRegressionComputer.setX(x);
        multipleLinearRegressionComputer.setxName(name);
        boolean finish = multipleLinearRegressionComputer.doMultipleLinearRegression();


        double predRSquared = multipleLinearRegressionComputer.getPredRSquared();
        System.out.println("predRSquared: "+predRSquared);
    }
}

class Parameter implements Comparable {
    public Parameter(String name, double tValue) {
        this.name = name;
        this.tValue = tValue;
    }

    String name;
    double tValue;

    public String getName() {
        return name;
    }

    public void setName(String name) {
        this.name = name;
    }

    public double gettValue() {
        return tValue;
    }

    public void settValue(double tValue) {
        this.tValue = tValue;
    }

    @Override
    public int compareTo(Object o) {
        return Double.compare(this.tValue, ((Parameter) o).tValue);
    }
}